Monte Carlo Methods: Lab 2

Point to note: an electronic copy of the Frenkel and Smit book is available through the library. This lab is based on case study 1 in chapter 3.4 of that book.

Lennard-Jones fluids

Use a Lennard-Jones potential inside a box size $[0,L]^3$ with a cut-off $r_c = L/2$,

$$ \begin{equation} U = \begin{cases} 4 \left[ \frac{1}{r^{12}} - \frac{1}{r^6} \right] & r < r_c \\ 0 & r > r_c. \end{cases} \end{equation} $$

Include tail corrections (that is, additional energy and pressure terms resulting from the particles outside the cutoff radius) as

$$ \begin{align} U^{\text{tail}} & = \frac{8 \pi \rho}{3} \left[ \frac{1}{3} \frac{1}{r_c^9} - \frac{1}{r_c^3} \right] \\ p^{\text{tail}} & = \frac{16 \pi \rho^2}{3} \left[ \frac{2}{3} \frac{1}{r_c^9} - \frac{1}{r_c^3} \right]. \end{align} $$

For each configuration we need to compute the pressure using

$$ \begin{equation} p = \frac{\rho}{\beta} + \frac{\text{Virial}}{V} \end{equation} $$

where

$$ \begin{equation} \text{Virial} = \sum_i \sum_{j > i} \vec{f}( \vec{r}_{ij} ) \cdot \vec{r}_{ij} \end{equation} $$

where, as usual, $\vec{r}_{ij}$ is the separation between the atoms, $\vec{r}_{ij} = \vec{r}_i - \vec{r}_j$, and the intermolecular force $\vec{f}$ is given by

$$ \begin{align} \vec{f}(\vec{r}_{ij}) &= - \nabla U \\ & = \begin{cases} 24 \left[ 2 \frac{1}{r^{14}} - \frac{1}{r^8} \right] \vec{r}_{ij} & r < r_c \\ \vec{0} & r > r_c \end{cases} \end{align} $$

Note that in the reduced coordinates $\beta = T^{-1}$.

Monte Carlo code

We will be using an $NTV$ approach, keeping the number of particles fixed ($N = 100$), the temperature fixed (either $T=2$ or $T=0.9$) and the volume fixed (indirectly, via the density $\rho = N / V = N L^{-3}$; use $\rho = a/10$ for $a = 1, \dots, 9$). You will need to take at least $10,000$ steps; $20,000$ is better, but in all cases you should test with a smaller number of particles and steps.

For reference we note the solutions, taken from Johnson, Zollweg and Gubbins for the pressures at $T=2$ are:

And for the pressures at $T=0.9$ are (noting that very little data is meaningful here - see the book):

Efficiency

Note that the sum over all particles scales as $n^2$ where $n$ is the number of particles. As the number of steps the algorithm will need to take will also scale as $n$, this makes the number of calculations at least as bad as $n^3$. This is expensive; if you try the naive approach then you'll have difficulty using more than 50 particles in a moderate time.

Instead we can note that, at each stage, the algorithm will move only one particle. Therefore, if we store not just the locations of the particles but also their pairwise separations, at each step we will only have to modify a small number of the separations. So we can store $r^2_{ij} = \vec{r}_{ij} \cdot \vec{r}_{ij}$ only, for $j > i$, and when perturbing particle $k$ we only need to update the separations $r^2_{ik}$ for $i<k$ and $r^2_{kj}$ for $k<j$.

This should significantly reduce the number of calculations done in each step.

In addition, note that for reasonable behaviour the acceptance rate should be $\sim 40\%$. This depends on the fractional perturbation distance $\Delta$; values $\sim 0.4$ are reasonable when $\rho \sim 0.1$, but values $\sim 0.02$ are reasonable when $\rho \sim 0.9$.

Results

• Check that the energy has converged to a "constant" state.
• Plot, for fixed temperature, $p$ vs $\rho$. The true solution is given on page 53 of the reference above.
• Plot a histogram of the energies to show that they follow the Boltzmann distribution.